1. Backgound

The models that are used here to derive the posterior samples were selected using rscripts loo comparison.


2. Set-up

Load all required packages for analysis.

library(tidyverse) # for data manipulation
library(pmthemes) # for ggplot themes
library(knitr) # for tables
library(sf)
## Warning: package 'sf' was built under R version 4.1.2
library(ggthemes)
library(viridis)
library(scales)
library(lubridate)
library(brms)
library(arsenal)
## Warning: package 'arsenal' was built under R version 4.1.2
library(knitr)
library(here)
## Warning: package 'here' was built under R version 4.1.1
library(spdep)
## Warning: package 'sp' was built under R version 4.1.2
library(broom)
library(RANN)
library(forcats)
library(cowplot)
library(bayesplot)
library(psych)
library(usethis)


3. Import datasets

Import data required for the analysis.

dat <- readRDS(here::here("data/dat.rds"))
dat <- sf::st_set_geometry(dat, NULL)

prev_model_rintc_31 <- readRDS(here::here("data/prev_model_rintc_31.rds"))
notif_model_rintc_25 <- readRDS(here::here("data/notif_model_rintc_25.rds"))


4. Check the prevalance model that was selected

summary(prev_model_rintc_31)
##  Family: zero_inflated_poisson 
##   Links: mu = log; zi = identity 
## Formula: n_prev_tbcases ~ 1 + scale_prop_adults_mean + (1 | cluster) + offset(log(tent_cxr_total)) 
##    Data: dat_scale (Number of observations: 72) 
## Samples: 3 chains, each with iter = 15000; warmup = 1000; thin = 1;
##          total post-warmup samples = 42000
## 
## Group-Level Effects: 
## ~cluster (Number of levels: 72) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.33      0.24     0.01     0.90 1.00    14078    16689
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                 -6.07      0.28    -6.63    -5.51 1.00    26216
## scale_prop_adults_mean    -0.06      0.08    -0.22     0.09 1.00    49468
##                        Tail_ESS
## Intercept                 28098
## scale_prop_adults_mean    28972
## 
## Family Specific Parameters: 
##    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## zi     0.18      0.13     0.01     0.46 1.00    27370    19772
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(prev_model_rintc_31)


5. Drawing the posterior distribution for the posterior model for prevalences

# Prevalence model

# extract the posterior samples from the prevalence model
# Remember need to use centred variables - describe clearly in the model code how you have centred.
post_tb_prev_m31 <- posterior_samples(prev_model_rintc_31)

# add the posterior intercepts to the random effects for each cluster

# mutate(across(contains("r_cluster"), ~((zi + (1-zi))*(exp(b_intercept + .x))))) %>%

prev_m_sum <- post_tb_prev_m31 %>%
  janitor::clean_names() %>%
  select(b_intercept, zi, contains("r_cluster")) %>%
  mutate(across(contains("r_cluster"), ~ (((1 - zi)) * (exp(b_intercept + .x))))) %>%
  select(-b_intercept, -zi) %>%
  rename_with(.fn = ~ stringr::str_replace_all(., "r_cluster", "p_prev")) %>%
  rename_with(.fn = ~ stringr::str_replace_all(., "_intercept", ""))


# have a look at the table
cluster_prevtotals <- select(dat, cluster, tent_cxr_total) %>%
  pivot_wider(names_from = cluster, values_from = tent_cxr_total) %>%
  rename_with(.fn = ~ stringr::str_replace_all(., "c", "p_prev_"))

# multiply by the cluster totals, to get the prevalent cases per cluster
# prev_m_sum <- prev_m_sum %>%
#                 mutate(across(everything(), ~.x * cluster_prevtotals[[cur_column()]]))
# prev_m_sum

prev_m_sum <- map2_df(prev_m_sum, cluster_prevtotals, ~ (.x * .y))

# posterior density plots for all clusters
# This gets the original prevalent cases per cluster from the study to compare with the model ones
cluster2019prev <- select(dat, cluster, n_prev_tbcases) %>%
  mutate(name = stringr::str_replace_all(cluster, "c", "p_prev_")) %>%
  mutate(name = fct_inorder(name))
prev_m_sum <- prev_m_sum %>%
  rename_all(.funs = funs(str_replace(., "c", "")))
## Warning: `funs()` was deprecated in dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
prev_m_sum <- prev_m_sum %>%
  select(num_range("p_prev_", range = 1:72))

prev_m_sum %>%
  pivot_longer(p_prev_1:p_prev_20) %>%
  select(name, value) %>%
  mutate(name = fct_inorder(name)) %>%
  ggplot() +
  geom_density(aes(x = value)) +
  geom_vline(
    data = cluster2019prev[1:20, ],
    aes(xintercept = n_prev_tbcases), color = "red"
  ) +
  facet_wrap(name ~ ., scales = "free") +
  theme_minimal()

prev_m_sum %>%
  pivot_longer(p_prev_1:p_prev_20) %>%
  select(name, value) %>%
  mutate(name = fct_inorder(name)) %>%
  ggplot() +
  geom_boxplot(aes(y = value, x = name)) +
  geom_hline(
    data = cluster2019prev[1:20, ],
    aes(yintercept = n_prev_tbcases),
    color = "red"
  ) +
  facet_wrap(name ~ ., scales = "free") +
  theme_minimal()

prev_m_sum %>%
  pivot_longer(p_prev_21:p_prev_41) %>%
  select(name, value) %>%
  mutate(name = fct_inorder(name)) %>%
  ggplot() +
  geom_density(aes(x = value)) +
  geom_vline(
    data = cluster2019prev[21:41, ], aes(xintercept = n_prev_tbcases),
    color = "red"
  ) +
  facet_wrap(name ~ ., scales = "free") +
  theme_minimal()

prev_m_sum %>%
  pivot_longer(p_prev_21:p_prev_41) %>%
  select(name, value) %>%
  mutate(name = fct_inorder(name)) %>%
  ggplot() +
  geom_boxplot(aes(y = value, x = name)) +
  geom_hline(
    data = cluster2019prev[21:41, ],
    aes(yintercept = n_prev_tbcases), color = "red"
  ) +
  facet_wrap(name ~ ., scales = "free") +
  theme_minimal()

prev_m_sum %>%
  pivot_longer(p_prev_42:p_prev_62) %>%
  select(name, value) %>%
  mutate(name = fct_inorder(name)) %>%
  ggplot() +
  geom_density(aes(x = value)) +
  geom_vline(
    data = cluster2019prev[42:62, ], aes(xintercept = n_prev_tbcases),
    color = "red"
  ) +
  facet_wrap(name ~ ., scales = "free") +
  theme_minimal()

prev_m_sum %>%
  pivot_longer(p_prev_42:p_prev_62) %>%
  select(name, value) %>%
  mutate(name = fct_inorder(name)) %>%
  ggplot() +
  geom_boxplot(aes(y = value, x = name)) +
  geom_hline(
    data = cluster2019prev[42:62, ],
    aes(yintercept = n_prev_tbcases), color = "red"
  ) +
  facet_wrap(name ~ ., scales = "free") +
  theme_minimal()

prev_m_sum %>%
  pivot_longer(p_prev_63:p_prev_72) %>%
  select(name, value) %>%
  mutate(name = fct_inorder(name)) %>%
  ggplot() +
  geom_density(aes(x = value)) +
  geom_vline(
    data = cluster2019prev[63:72, ], aes(xintercept = n_prev_tbcases),
    color = "red"
  ) +
  facet_wrap(name ~ ., scales = "free") +
  theme_minimal()

prev_m_sum %>%
  pivot_longer(p_prev_63:p_prev_72) %>%
  select(name, value) %>%
  mutate(name = fct_inorder(name)) %>%
  ggplot() +
  geom_boxplot(aes(y = value, x = name)) +
  geom_hline(
    data = cluster2019prev[63:72, ],
    aes(yintercept = n_prev_tbcases), color = "red"
  ) +
  facet_wrap(name ~ ., scales = "free") +
  theme_minimal()


6. Drawing the posterior distribution for the posterior model for notifications

summary(notif_model_rintc_25)
##  Family: poisson 
##   Links: mu = log 
## Formula: total_confirmed ~ 1 + scale_prop_adults_mean + scale_perc_never_primary_mean + scale_clinic_distance_1km + tb_year + (1 | cluster) + offset(log(population)) 
##    Data: dat_scale_ln (Number of observations: 360) 
## Samples: 3 chains, each with iter = 15000; warmup = 1000; thin = 1;
##          total post-warmup samples = 42000
## 
## Group-Level Effects: 
## ~cluster (Number of levels: 72) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.31      0.04     0.24     0.39 1.00    15832    25699
## 
## Population-Level Effects: 
##                               Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                        -7.58      0.08    -7.75    -7.42 1.00
## scale_prop_adults_mean           -0.04      0.02    -0.07     0.00 1.00
## scale_perc_never_primary_mean    -0.02      0.01    -0.04    -0.01 1.00
## scale_clinic_distance_1km        -0.25      0.06    -0.37    -0.13 1.00
## tb_year2015                       1.06      0.08     0.91     1.21 1.00
## tb_year2016                       1.07      0.08     0.92     1.22 1.00
## tb_year2017                       0.92      0.08     0.77     1.07 1.00
## tb_year2018                       0.20      0.09     0.03     0.37 1.00
##                               Bulk_ESS Tail_ESS
## Intercept                        18133    25410
## scale_prop_adults_mean           16997    24775
## scale_perc_never_primary_mean    18218    26704
## scale_clinic_distance_1km        17552    25411
## tb_year2015                      25266    30939
## tb_year2016                      25594    30585
## tb_year2017                      25767    30568
## tb_year2018                      28496    31935
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(notif_model_rintc_25)


6. Drawing the posterior distribution for the posterior model for notifications

# now do the same for the notification data
post_tb_notifs_m25 <- posterior_samples(notif_model_rintc_25)

# add the posterior intercepts to the random effects for each cluster
notif_m_sum <- post_tb_notifs_m25 %>%
  janitor::clean_names() %>%
  select(b_intercept, contains("r_cluster")) %>%
  mutate(across(contains("r_cluster"), ~ .x + b_intercept)) %>%
  select(-b_intercept) %>%
  rename_with(.fn = ~ stringr::str_replace_all(., "r_cluster", "p_notif")) %>%
  rename_with(.fn = ~ stringr::str_replace_all(., "_intercept", "")) %>%
  mutate(across(everything(), ~ exp(.x)))

# have a look at the table
# Get the cluster totals
cluster_totals <- select(dat, cluster, total) %>%
  pivot_wider(names_from = cluster, values_from = total) %>%
  rename_with(.fn = ~ stringr::str_replace_all(., "c", "p_notif_"))
# convert to notifications per cluster
# notif_m_sum <- notif_m_sum %>%
#                   mutate(across(everything(), ~.x * cluster_totals[[cur_column()]]))

notif_m_sum <- map2_df(notif_m_sum, cluster_totals, ~ (.x * .y))
notif_m_sum
# posterior density plots for all clusters
# Get the number of TB notifications for 2019 to compare with the model one's
cluster2019notifs <- select(dat, cluster, total_tbcases_2019) %>%
  mutate(name = stringr::str_replace_all(cluster, "c", "p_notif_")) %>%
  mutate(name = fct_inorder(name))


notif_m_sum <- notif_m_sum %>%
  rename_all(.funs = funs(str_replace(., "c", "")))

notif_m_sum <- notif_m_sum %>%
  select(num_range("p_notif_", range = 1:72))


notif_m_sum %>%
  pivot_longer(p_notif_1:p_notif_20) %>%
  select(name, value) %>%
  mutate(name = fct_inorder(name)) %>%
  ggplot() +
  geom_density(aes(x = value)) +
  geom_vline(
    data = cluster2019notifs[1:20, ],
    aes(xintercept = total_tbcases_2019), color = "red"
  ) +
  facet_wrap(name ~ ., scales = "free") +
  theme_minimal()

notif_m_sum %>%
  pivot_longer(p_notif_1:p_notif_20) %>%
  select(name, value) %>%
  mutate(name = fct_inorder(name)) %>%
  ggplot() +
  geom_boxplot(aes(y = value, x = name)) +
  geom_hline(
    data = cluster2019notifs[1:20, ],
    aes(yintercept = total_tbcases_2019), color = "red"
  ) +
  facet_wrap(name ~ ., scales = "free") +
  theme_minimal()

notif_m_sum %>%
  pivot_longer(p_notif_21:p_notif_41) %>%
  select(name, value) %>%
  mutate(name = fct_inorder(name)) %>%
  ggplot() +
  geom_density(aes(x = value)) +
  geom_vline(
    data = cluster2019notifs[21:41, ], aes(xintercept = total_tbcases_2019),
    color = "red"
  ) +
  facet_wrap(name ~ ., scales = "free") +
  theme_minimal()

notif_m_sum %>%
  pivot_longer(p_notif_21:p_notif_41) %>%
  select(name, value) %>%
  mutate(name = fct_inorder(name)) %>%
  ggplot() +
  geom_boxplot(aes(y = value, x = name)) +
  geom_hline(
    data = cluster2019notifs[21:41, ],
    aes(yintercept = total_tbcases_2019), color = "red"
  ) +
  facet_wrap(name ~ ., scales = "free") +
  theme_minimal()

notif_m_sum %>%
  pivot_longer(p_notif_42:p_notif_62) %>%
  select(name, value) %>%
  mutate(name = fct_inorder(name)) %>%
  ggplot() +
  geom_density(aes(x = value)) +
  geom_vline(
    data = cluster2019notifs[42:62, ], aes(xintercept = total_tbcases_2019),
    color = "red"
  ) +
  facet_wrap(name ~ ., scales = "free") +
  theme_minimal()

notif_m_sum %>%
  pivot_longer(p_notif_42:p_notif_62) %>%
  select(name, value) %>%
  mutate(name = fct_inorder(name)) %>%
  ggplot() +
  geom_boxplot(aes(y = value, x = name)) +
  geom_hline(
    data = cluster2019notifs[42:62, ],
    aes(yintercept = total_tbcases_2019), color = "red"
  ) +
  facet_wrap(name ~ ., scales = "free") +
  theme_minimal()

notif_m_sum %>%
  pivot_longer(p_notif_63:p_notif_72) %>%
  select(name, value) %>%
  mutate(name = fct_inorder(name)) %>%
  ggplot() +
  geom_density(aes(x = value)) +
  geom_vline(
    data = cluster2019notifs[63:72, ], aes(xintercept = total_tbcases_2019),
    color = "red"
  ) +
  facet_wrap(name ~ ., scales = "free") +
  theme_minimal()

notif_m_sum %>%
  pivot_longer(p_notif_63:p_notif_72) %>%
  select(name, value) %>%
  mutate(name = fct_inorder(name)) %>%
  ggplot() +
  geom_boxplot(aes(y = value, x = name)) +
  geom_hline(
    data = cluster2019notifs[63:72, ],
    aes(yintercept = total_tbcases_2019), color = "red"
  ) +
  facet_wrap(name ~ ., scales = "free") +
  theme_minimal()


7. Prevalence to notification ratios

# Now to derive the prevalance to notification ratios
# First step is to scale the posterior draws to preve and not per 250000 population
# as in p1 and n1

p1 <- prev_m_sum %>%
  mutate(across(everything(), ~ ((.x / cluster_prevtotals[[cur_column()]]) * 100000))) %>%
  pivot_longer(p_prev_1:p_prev_72) %>%
  mutate(cluster = parse_number(name)) %>%
  rename(p_prev = value) %>%
  select(-name)

n1 <- notif_m_sum %>%
  mutate(across(everything(), ~ ((.x / cluster_totals[[cur_column()]]) * 100000))) %>%
  pivot_longer(p_notif_1:p_notif_72) %>%
  # mutate(cluster = parse_number(name)) %>%
  rename(p_notif = value) %>%
  select(-name)


# join together the posterior predictive draws for the prev model and notifs model.
p_prev31_notif25_rintc <- bind_cols(p1, n1)

# Take the ratio of the prevalence and notifications,
p_prev31_notif25_rintc <- p_prev31_notif25_rintc %>%
  mutate(prev_to_notif_ratio = p_prev / p_notif)

# look at this table
tidybayes::mean_qi(p_prev31_notif25_rintc$prev_to_notif_ratio)
p_prev31_notif25_rintc
# now group by cluster and summarise, here using medians and 50% intervals (but you can change this)
library(tidybayes)
## 
## Attaching package: 'tidybayes'
## The following objects are masked from 'package:brms':
## 
##     dstudent_t, pstudent_t, qstudent_t, rstudent_t
cluster_ests_50_prev31_notif25_rintc <- p_prev31_notif25_rintc %>%
  group_by(cluster) %>%
  mean_qi(prev_to_notif_ratio, .width = 0.50)

# look at this table
cluster_ests_50_prev31_notif25_rintc
# plot the prevalence to notification ratios per cluster
cluster_ests_50_prev31_notif25_rintc %>%
  ggplot() +
  geom_pointrange(aes(x = cluster, y = prev_to_notif_ratio, ymin = .lower, ymax = .upper)) +
  geom_hline(aes(yintercept = 1), colour = "red", linetype = "dashed") +
  labs(y = "Prevalence to notification ratio (50% CI)")

# 95% CI
cluster_ests_95_prev31_notif25_rintc <- p_prev31_notif25_rintc %>%
  group_by(cluster) %>%
  mean_qi(prev_to_notif_ratio, .width = 0.95)

# The overall 95% estimates
p_prev31_notif25_rintc %>%
  mean_qi(prev_to_notif_ratio, .width = 0.95)
# look at this table
cluster_ests_95_prev31_notif25_rintc
# plot the prevalence to notification ratios per cluster
cluster_ests_95_prev31_notif25_rintc %>%
  ggplot() +
  geom_pointrange(aes(x = cluster, y = prev_to_notif_ratio, ymin = .lower, ymax = .upper)) +
  geom_hline(aes(yintercept = 1), colour = "red", linetype = "dashed") +
  labs(y = "Prevalence to notification ratio (95% CI)")

saveRDS(object = post_tb_prev_m31, file = here::here("data/post_tb_prev_m31.rds"))
saveRDS(object = post_tb_notifs_m25, file = here::here("data/post_tb_notifs_m25.rds"))

saveRDS(object = p_prev31_notif25_rintc, file = here::here("data/p_prev31_notif25_rintc.rds"))
saveRDS(object = cluster_ests_50_prev31_notif25_rintc, file = here::here("data/cluster_ests_50_prev31_notif25_rintc.rds"))
saveRDS(object = cluster_ests_95_prev31_notif25_rintc, file = here::here("data/cluster_ests_95_prev31_notif25_rintc.rds"))


8. Reproducibility

This reproduction of the analysis was run by:

keyName value
sysname Windows
release 10 x64
version build 19042
nodename LP256
machine x86-64
login mkhundi
user mkhundi
effective_user mkhundi

Analysis was run at 2022-02-03 16:04:17, and using the following Session Info:

R version 4.1.0 (2021-05-18)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.1252 
[2] LC_CTYPE=English_United Kingdom.1252   
[3] LC_MONETARY=English_United Kingdom.1252
[4] LC_NUMERIC=C                           
[5] LC_TIME=English_United Kingdom.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] tidybayes_2.3.1   usethis_2.0.1     psych_2.1.6       bayesplot_1.8.1  
 [5] cowplot_1.1.1     RANN_2.6.1        broom_0.7.8       spdep_1.1-8      
 [9] spData_0.3.10     sp_1.4-6          here_1.0.1        arsenal_3.6.3    
[13] brms_2.15.0       Rcpp_1.0.7        lubridate_1.7.10  scales_1.1.1     
[17] viridis_0.6.1     viridisLite_0.4.0 ggthemes_4.2.4    sf_1.0-5         
[21] knitr_1.33        pmthemes_0.2.0    forcats_0.5.1     stringr_1.4.0    
[25] dplyr_1.0.7       purrr_0.3.4       readr_1.4.0       tidyr_1.1.3      
[29] tibble_3.1.2      ggplot2_3.3.5     tidyverse_1.3.1  

loaded via a namespace (and not attached):
  [1] readxl_1.3.1         backports_1.2.1      plyr_1.8.6          
  [4] igraph_1.2.6         svUnit_1.0.6         splines_4.1.0       
  [7] crosstalk_1.1.1      rstantools_2.1.1     inline_0.3.19       
 [10] digest_0.6.27        htmltools_0.5.1.1    rsconnect_0.8.18    
 [13] gdata_2.18.0         fansi_0.5.0          magrittr_2.0.1      
 [16] modelr_0.1.8         gmodels_2.18.1       RcppParallel_5.1.4  
 [19] matrixStats_0.59.0   xts_0.12.1           prettyunits_1.1.1   
 [22] colorspace_2.0-2     rvest_1.0.0          ggdist_2.4.1        
 [25] haven_2.4.1          xfun_0.24            callr_3.7.0         
 [28] crayon_1.4.2         jsonlite_1.7.2       lme4_1.1-27.1       
 [31] zoo_1.8-9            glue_1.4.2           gtable_0.3.0        
 [34] emmeans_1.6.1        V8_4.0.0             distributional_0.2.2
 [37] pkgbuild_1.2.0       rstan_2.21.2         abind_1.4-5         
 [40] mvtnorm_1.1-2        DBI_1.1.2            miniUI_0.1.1.1      
 [43] xtable_1.8-4         tmvnsim_1.0-2        units_0.7-2         
 [46] proxy_0.4-26         stats4_4.1.0         StanHeaders_2.21.0-7
 [49] DT_0.18              htmlwidgets_1.5.3    httr_1.4.2          
 [52] threejs_0.3.3        arrayhelpers_1.1-0   ellipsis_0.3.2      
 [55] farver_2.1.0         pkgconfig_2.0.3      loo_2.4.1           
 [58] deldir_0.2-10        sass_0.4.0           dbplyr_2.1.1        
 [61] janitor_2.1.0        utf8_1.2.1           labeling_0.4.2      
 [64] tidyselect_1.1.1     rlang_0.4.11         reshape2_1.4.4      
 [67] later_1.2.0          munsell_0.5.0        cellranger_1.1.0    
 [70] tools_4.1.0          cli_3.0.0            generics_0.1.0      
 [73] ggridges_0.5.3       evaluate_0.14        fastmap_1.1.0       
 [76] yaml_2.2.1           processx_3.5.2       fs_1.5.0            
 [79] nlme_3.1-152         mime_0.12            projpred_2.0.2      
 [82] xml2_1.3.2           compiler_4.1.0       shinythemes_1.2.0   
 [85] rstudioapi_0.13      gamm4_0.2-6          curl_4.3.2          
 [88] e1071_1.7-9          reprex_2.0.0         bslib_0.2.5.1       
 [91] stringi_1.6.1        highr_0.9            ps_1.6.0            
 [94] Brobdingnag_1.2-6    lattice_0.20-44      Matrix_1.3-3        
 [97] classInt_0.4-3       nloptr_1.2.2.2       markdown_1.1        
[100] shinyjs_2.0.0        vctrs_0.3.8          LearnBayes_2.15.1   
[103] pillar_1.6.4         lifecycle_1.0.1      jquerylib_0.1.4     
[106] bridgesampling_1.1-2 estimability_1.3     raster_3.5-11       
[109] httpuv_1.6.1         R6_2.5.1             promises_1.2.0.1    
[112] KernSmooth_2.23-20   gridExtra_2.3        codetools_0.2-18    
[115] boot_1.3-28          colourpicker_1.1.0   MASS_7.3-54         
[118] gtools_3.9.2         assertthat_0.2.1     rprojroot_2.0.2     
[121] withr_2.4.3          mnormt_2.0.2         shinystan_2.5.0     
[124] expm_0.999-6         mgcv_1.8-35          parallel_4.1.0      
[127] hms_1.1.1            terra_1.4-22         grid_4.1.0          
[130] coda_0.19-4          class_7.3-19         minqa_1.2.4         
[133] snakecase_0.11.0     rmarkdown_2.9        shiny_1.6.0         
[136] base64enc_0.1-3      dygraphs_1.1.1.6